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Abstract 

A three-dimensional (3D) lattice Boltzmann method based on central moments is 
derived. Two main elements are the local attractors in the collision term and the 
source terms representing the effect of external and/or self-consistent internal forces. 
For suitable choices of the orthogonal moment basis for the three-dimensional, 
twenty seven velocity (D3Q27), and, its subset, fifteen velocity (D3Q15) lattice 
models, attractors are expressed in terms of factorization of lower order moments 
as suggested in an earlier work; the corresponding source terms are specified to 
correctly influence lower order hydrodynamic fields, while avoiding aliasing effects 
for higher order moments. These are achieved by successively matching the cor- 
responding continuous and discrete central moments at various orders, with the 
final expressions written in terms of raw moments via a transformation based on 
the binomial theorem. Furthermore, to alleviate the discrete effects with the source 
terms, they are treated to be temporally semi-implicit and second-order, with the 
implicitness subsequently removed by means of a transformation. As a result, the 
approach is frame-invariant by construction and its emergent dynamics describing 
fully 3D fluid motion in the presence of force fields is Galilean invariant. Numerical 
experiments for a set of benchmark problems demonstrate its accuracy. 
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1 Introduction 



The use of discrete velocity models based on kinetic theory is a powerful 
theoretical approach and forms the basis of a modern computational method 
for fluid mechanics. While the work of Broadwell [1] represents an early effort 
in this direction, careful exploitation of symmetries and local conservation 
laws to construct such models for discrete configuration spaces underpinned 
the recent approaches, starting from the work of Frisch et al [2J. The latter 
led to the development of the lattice Boltzmann method (LBM) [3], albeit 
without any direct connection to kinetic theory in its initial stages. Indeed, 
formal demonstration of this approach as a simplified model for the continuous 
Boltzmann equation [UEJEj, provided much impetus for recent developments, 
particularly for complex fluids [7,8. 9j and for representation beyond continuum 
description [10] . among others (see [TTlll2|ll3|fH] for general reviews on the 
LBM). 

The basic procedure involved in the LBM is represented by the synchronous 
free-streaming of particle distribution functions along discrete directions fol- 
lowed by collision, represented as a relaxation process. The latter has major 
influence on the physical fidelity as well as numerical stability. A popular 
approach is based on the single-relaxation time (SRT) model [T5][T6] . While 
it is successful in many applications, it is prone to numerical instability for 
situations with relatively low viscosities and is inadequate for representing 
certain physical phenomena (e.g. viscoelasticity and thermal transport) and 
in correctly accounting for kinetic layers near boundaries. In contrast, the use 
of multiple relaxation time (MRT) models [T7], which are simplified versions 
of the relaxation LBM [T8|IT9] . have addressed these aspects significantly. Its 
characteristic feature is that the relaxation process is carried out in moment 
space [20]. In particular, the relaxation times for the kinetic modes can be 
independently adjusted by means of a linear stability analysis to improve nu- 
merical stability [2~Tf22] . Furthermore, based on the notion of duality between 
hydrodynamic and kinetic modes, a procedure for construction of matrix based 
LBM has been proposed recently [23] • From a different standpoint, non-linear 
stability can be enforced with a discrete H-theorem locally in the collision 
step using the SRT model [2^|25II26[ . In this Entropic LBM, minimization of 
a convex H-function with hydrodynamic conservation constraints yields tran- 
scendental local attractors. It was also shown that the choice of the H-function 
in this context can be determined by enforcing Galilean invariance [27]. The 
construction based on the minimization of a convex function has been general- 
ized to include a larger set of constraints that includes second-order moments 
yielding quasi-equilibrium attractors and thus allowing for a two-relaxation 
time Entropic LBM via a continuous H-theorem |28|I29] . A theoretical ba- 
sis for such an approach based on factorization symmetry considerations has 
been presented in [30] • This work, along with others [3~T|32f33j . also provides 
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rational procedures for constructing higher order models. 

For general applicability of models and numerical schemes, it is necessary that 
their description of fluid behavior be the same in all inertial frames of reference. 
This important physical requirement of Galilean invariance, when not met can 
also lead to numerical instability in the context of the LBM. The latter arises 
from the fact that the degeneracies due to the fmiteness of the standard lattice 
velocity sets can lead to linear dependence of higher order moments in terms 
of those at lower orders, which, in turn, can result in negative dependence of 
viscosity on fluid velocity [HI]. Thus, it becomes necessary to consider large 
lattice velocity sets, which, however, by themselves do not guarantee in strictly 
observing Galilean invariance, as they can only lead to kinematically complete 
models [35] • Proper selection of the collision model provides the sufficient 
or the dynamically complete condition in this regard to recover the correct 
physics, such as the Navier-Stokes equations. This can be seen by the use of 
unwieldy fitting of parameters [34] or elaborate construction procedures [32] 
for the attractors in the collision model with such extended lattice sets. Thus, 
the collision process still needs to be carefully designed and has an important 
role to play in the proper observation of Galilean invariance. 

In this context, relaxation in a moving frame of reference, i.e. in terms of mo- 
ments obtained by shifting the particle velocity with the local hydrodynamic 
velocity, or central moments provides a natural setting and a simple con- 
struction procedure to maintain Galilean invariance for a given velocity set. 
That is, the relaxation process is constructed to observe inertial frame invari- 
ance to a degree as permitted by the chosen lattice velocity set. We consider 
this specific meaning when we use the term Galilean invariance in this pa- 
per. Also, when different central moments are relaxed at different rates, i.e. 
formulated as a MRT model, it can enhance numerical stability by provid- 
ing additional numerical dissipation similar to standard MRT models based 
on raw moments. It is noted that the ideas and procedures based on central 
moment relaxation are not restricted to standard lattice velocity sets, but can 
be used for any extended or kinematically complete velocity sets as well. The 
central moment approach exhibits a cascaded structure, which was shown to 
be equivalent to considering a generalized equilibrium in the lattice or rest 
frame of reference [3B]- Furthermore, to further improve the physical fidelity, 
the local attractors for the central moments given in terms of their factor- 
ization into lower order moments has been proposed [39J. To incorporate the 
effect of force fields, which are important for numerous physical applications, 
a new approach for the source terms based on central moments was recently 
developed for a two-dimensional (2D) lattice [ID] . In addition, a detailed the- 
oretical basis for the central moment method, including a consistency analysis 
of the emergent fluid motion, was also provided [40] . 

The objective of this work is the derivation and validation of a 3D central mo- 
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ment lattice Boltzmann method, with a particular focus on deriving Galilean 
invariant source terms, which are important, for example, in situations involv- 
ing multiphase/multicomponent flows or turbulence modeling. In this regard, 
three-dimensional, twenty seven velocity (D3Q27) and its minimal subset, i.e. 
fifteen velocity (D3Q15) velocity lattices that can recover Navier-Stokes be- 
havior are considered, and the overall procedure and notations used in [ID] are 
adopted. The D3Q27 lattice is chosen so that our results provide the forcing 
scheme based on central moments to the overall formulation considered in [36] . 
It is noted that the notations and the details provided in that work [36] are 
cumbersome even for the collision model without forcing for implementation. 
On the other hand, in practice, the computational complexity is considerably 
reduced with the use of the D3Q15 lattice. Hence, the details with both the 
lattices are provided, with the smaller lattice set used in most of the compu- 
tations in our validation studies. The overall procedure is as follows. Starting 
from suitable choices of the orthogonal moment basis for these lattice veloc- 
ity sets, the continuous and discrete central moments of the local attractors 
and source terms at different orders are successively matched. The results are 
then transformed in terms of raw moments by means of the binomial theo- 
rem. To maintain physical coherence for the discrete velocity set, factorized 
local attractors for higher order central moments and temporally second-order 
accurate treatment of source terms are considered. This construction yields 
Galilean invariant representation of 3D fluid dynamics in the presence of gen- 
eral external or self-consistent internal forces. The computational approach 
thus derived is then assessed by comparison of its results for a set of canonical 
problems involving forcing for which analytical solutions are available. 



The paper is organized as follows. Its main body containing the derivation 
focuses only on the essential steps involved in the derivation, choosing the 
D3Q27 lattice as an example, with the attendant details presented in various 
appendices (see Appendices [A 



lattice is presented in Appendix 



G 



the computational scheme for the D3Q15 
G). Section [2] discusses the choice made for the 
orthogonal moment basis corresponding to the D3Q27 lattice. The ansatz for 
the continuous central moments for the distribution functions, local attractors 
and sources due to the force fields are presented in Sees. [3] Section [4] provides 
the corresponding 3D lattice Boltzmann equation (LBE) with source terms 
based on central moments. Various discrete central moments needed for the 
construction of the central moment method are defined and the matching 
principle to preserve Galilean invariance is stated in Sec. [5j Section [6] obtains 
expressions for various discrete raw moments using the matching principle via 
the binomial theorem, including the derivation of the source terms in particle 
velocity space. The construction of the collision kernel is presented in Sec. [7] 
and the overall procedure of the central moment LBM is provided in Sec. [8] 
Validation studies involving various canonical problems are discussed in Sec. [9] 



The conclusions are finally summarized in Sec. 10 
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2 Selection of Moment Basis 



We now discuss the moment basis, which is an important element on which the 
central moment LBM is constructed, corresponding to the three-dimensional, 
twenty seven velocity (D3Q27) lattice model (see Fig. [I]). The particle velocity 
for this lattice model ~~i ' a is given by 



wv 



(0,0,0), a = 

(±1,0,0),(0,±1,0),(0,0,±1), a = l,-..,6 
(±1, ±1, 0), (±1, 0, ±1), (0, ±1, ±1), a = 7, • • • , 18 
(±1,±1,±1), a = 19,---,26 



(1) 




Fig. 1. Three-dimensional, twenty seven particle velocity (D3Q27) lattice. 

For convenience, as in [40J, we use Dirac's bra-ket notion to represent the 
basis vectors, and Greek and Latin subscripts for particle velocity directions 
and Cartesian coordinate directions, respectively. By definition, the moments 
in the LBM are discrete integral properties of the distribution function f a , i.e. 
Y?a=o e ax e ay e azfa, where m, u and p are integers, in 3D. As a result, we begin 
with the following twenty-seven non-orthogonal independent basis vectors ob- 
tained by combining monomials e^e™ e p az and arranged in an ascending order. 
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First, the nominal basis for the conserved (zeroth and first order) moments 
follows immediately: 



To) = 


\P) = 


Ti) = 




Ti) = 




T 3 ) = 





The basis for second-order moments are chosen such that it allows correct rep- 
resentation of the momentum flux (based on the Maxwell distribution, see be- 
low in Sec. |3| in the hydrodynamic equations. Three off-diagonal components 
(|T 4 )-|T 6 )) and three diagonal components (\T 7 )-\T 9 )) with max(m,n,p) = 1 
and max(m, n,p) = 2, respectively, while satisfying m + n + p = 2 are consid- 
ered: 



\T*) = 


\&ax&ay) > 


\T 5 ) = 




\n) = 




\T 7 ) = 


1 2 2 \ 
IP — P ) 
l^ax ay/ ' 


\Ts) = 


1 2 2 \ 
\ e ax e az) ' 


\T 9 ) = 


1 2 _|_ 2 _|_ 2 \ 
\ e ax ' e ay ' e az) 



The following six third-order basis vectors for moments are chosen (|T 10 )- 
IT15) with max(m,n,p) = 2 and |Tie) with max(m,n,p) = 1, while satisfying 
m + n + p = 3): 



Ti ) = 


\ e ax e ay + e «x e Q2 ) 5 


Tu) = 




T 12 ) = 


l e «:c e az ~l~ &ay^az) j 


T 13 ) = 


1 2 2 \ 


T u ) = 


1 2 2 \ 
l e aa; e aj/ e &y e az) ) 


T 15 ) = 


1 2 2 \ 
l e Qx e az e ay e az) j 


Ti 6 ) = 


| &ax&ay&otz) ■ 



For the fourth-order basis vectors, we consider the following six of them (IT17)- 
IT22) with max(m, n,p) = 2 for m + n + p = 4) : 



1^7) 
l^is) 

I^19> 



I 2 2 

P P 
I ^ax^ay 

I 2 2 

P P 

I 2 2 

P P 



1 2 2 
4- p p 

1 ^ax^az 



+ e 2 e 2 ) 



,2 



2 

az 



,2 



2 ) 



2 2 \ 
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|^2o) — \Cax e ay e az) , 
|-^~2l) I &ax& a y^otz) , 
1^22} — \^ax e ay e az) ■ 

Finally, three fifth-order basis vectors (IT23)— IT25)) and one sixth-order basis 
vector (|T2e)) are considered to complete moment basis corresponding to the 
D3Q27 model. In the above, in each case max(m, n,p) = 2, with m + n+p = 5 
and m + n + p = 6, respectively. Thus, 



1^23) 

l^24> 
1^25) 
1^26) 



I 2 2 \ 

\ e ax^ a y e azl 

I 2 2 
I 2 2 

e~„e~„,e 



>>x^ay^az) J 



I 2 2 2 



(2) 



Note that due to the finiteness of the particle velocity set, higher order lon- 
gitudinal moments, i.e. |e^) with m > 3 are eliminated from consideration in 
the above. The components of the basis vectors for the conserved moments 
corresponding to Eq. may be written as 

\p) = ||1^| > = (1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1,1,1,1,1,1)+, 

\e ax ) = (0, 1, -1, 0, 0, 0, 0, 1, -1, 1, -1, 1, -1, 1, -1, 0, 0, 0, 0, 1, -1, 

1,-1,1,-1,1,-1)+, 
\e ay ) = (0, 0, 0, 1, -1, 0, 0, 1, 1, -1, -1, 0, 0, 0, 0, 1, -1, 1, -1, 1, 1, 

-1,-1,1,1,-1,-1)+, 
\e az ) = (0,0, 0,0, 0,1, -1,0, 0,0, 0,1, 1,-1, -1,1, 1,-1, -1,1,1, 
1,1,-1,-1,-1,-1)+. 

Here and henceforth, the superscript '+' represents the transpose operator. The 
next key step is to transform the above non-orthogonal nominal basis set into 
an equivalent orthogonal set that would allow an efficient implementation [T7] . 
This is accomplished by means of the standard Gram-Schmidt procedure for 
the above arrangement, i.e. in the increasing order of the monomials of the 
products of the Cartesian components of the particle velocities. As a result 
the components of the orthogonal basis vectors are given by 



K ) = 


\p) = Wta\°), 


Ki) = 


\&ax) , 


K 2 ) = 




K 3 ) = 
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\Ka) 


\&ax&ay) 1 






\&ax&az) 1 






\&ay^az) 1 




1^7) 


1 2 2 \ 
\^ax ^ayl ) 




l^8> 


= | e 2 + e 2 +e 2 )- 

\^ax 1 ay 1 ^azl 


3 \ e az) > 


1^9) 


1 2 , 2 , 2 \ 
= P -4- P -4- P > — 

l^oa; 1 ay 1 ^azl 


2|p), 


#10) 


_ q 1 2 1 2 \ 


4 | Pax) , 


K n ) 


q 1 2 1 2 \ 

— " \ e ax e ay ~+~ e ay e az ) 


4 |6 a y) , 


K l2 ) 


q 1 2 _|_ 2 \ 

~~ " I e aa; e a2 e ay e otz) 


4 |pq. 2 ) , 


K 13 ) 


- 1 2 2 \ 




K u ) 


— 1 2 2 \ 
— \ e ax e ay e ay e a z) ' 




K 15 ) 


_ 1 2 2 \ 
~~ \ e ax e az e ay e az) , 




K 16 ) 






K17) 


1 2 2 1 2 2 1 2 2 \ 
— 1 P P -1- P P -1- P P ) 


K 18 ) 


1 2 2 1 2 2 
— " l e ox e aj/ ~r e oi e az 


- 2p 2 p 2 ) - 

^ay^azl 


K 19 ) 


1 2 2 2 2 \ 
— 1 P P P P } 

rac aj/ ^ax^azl 


' 2 \e 2 ay - « 


K20) 


3 |6 Q , a .6 Q ,^P Q2 ) 2 |p 


ay&az) 5 


K 2 i) 


3 |6 q , 2: p q ,^p Q2 ) 2 |p 


ax&az) 1 


K22) 


3 |6 Q , :r 6ay6 Q2 ) 2 |p 


ax&ay) 1 


K23) 


= 9 |6 Qa: 6 Q , ?/ P Q2 ) — 6 |p 


otx&ay 6ax 


K24) 


= 9 6 Qa .p Q ,j / p Q2 ) — 6 |e 


ax e ay + &ay 


K 25 ) 


~~ 9 | e cK2! e a2/ e a:z) 6 | e 


2 1 2 


K 26 ) 


— 27 Ip 2 p 2 p 2 \ — 18 

* ' \^ax^ay^azl ^ 


1 2 2 _|_ 

V^ax^ay ' e 



,2 
'ay 



2|2p 



2 

(i:c 



cm/ 



2 2 \ 

3 P ; 
"ay^az/ 



+ 12|eL + e 2 +eL)-8|p). 



(3) 



This can be explicitly written in terms of a orthogonal matrix of moment basis 
K given by 



K = [\K ) , IK,) , \K 2 ) , \K 3 ) , \K 4 ) , \K 5 ) , \K 6 ) , \K 7 ) , \K 8 ) 

\K Q ) , \K 10 ) , \K n ) , \K 12 ) , \K 13 ) , \K U ) , \K 15 ) , \K 16 ) , \K l7 }} . 

\K 18 ) , \K W ) , \K 20 ) , \K 21 ) , \K 22 ) , \K 23 ) , \K 2A ) , \K 25 ) , \K 2& )\ (4) 

whose components are presented in Appendix |XJ Note that unlike the stan- 
dard MRT formulation based on raw moments [22], which orders the basis 
vectors by considering the character of moments, i.e. increasing powers of 
their tensorial orders (i.e. scalars, vectors, tensors of different ranks,. . . ), the 
central moment basis vectors considered here are ordered according to their 
ascending powers (i.e. zeroth order moment, first order moments, second order 
moments,. . . ). Furthermore, the details of the basis vectors considered in this 
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paper are different from those provided in [36] 



3 Continuous Central Moments: Distribution Function, its Local 
Attractor and Forcing 



The central moment LBM, which is defined at the discrete level, should pre- 
serve certain continuous integral properties of the distribution function / given 
in terms of its central moments, i.e. those shifted by the macroscopic fluid ve- 
locity. In this regard, we first define continuous central moment of / of order 
(m + n + p) as 

oo oo oo 

U xmynzP = J J J - u x ) m (i v - u y ) n (£ z - u z ) p d£ x di y d£ z . (5) 



Here, and in the rest of this paper, the use of "hat" over a symbol represents 
quantities in the space of moments. The effect of collision is to relax the 
distribution function, or equivalently, its central moments, towards its local 
attractor. The corresponding central moment local attractor may be written 

as 

oo oo oo 

/ / / f at ^- u -) m ^y- u y) n ^- u ^ d ^y d ^- ( 6 ) 



n 



at 

m y U gp 



Here f at is as yet unknown, and its effect on the dynamics will be determined 
in what follows. Similarly, the continuous central moments due to sources may 
be written as 



oo oo oo 



?l mynzP = J J J Af F (£ x - u x ) m (£ y - u y ) n (i z - u z ) p d£ x d£ y d£ z , (7) 



— OO — oo — oo 



where Af F is the change in the distribution function due to forcing, which 
will be specified later. One possibility is to consider the local Maxwellian as 
the attractor [36]. That is, consider 

f M ee I" (/', Tt, („ ?„. f J = A exp 



2c? 
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where c 2 s = 1/3, which yields corresponding continuous Maxwellian central 
moments as 



oo oo oo 



ft£w=/ / j f M {^-n x ) m {i y -u y ) n ^ z -u z ydi x d^di z . (9) 



-oo — oo — oo 



By virtue of the fact that f M being an even function, Yl^i, y „ zP ^ when 
m, n and p are even and H%L ynzP = when m or n or p is odd. Here and 
henceforth, the subscripts x m y n z p mean xxx ■ ■ ■ TO-times, yyy ■ ■ ■ n-times and 
zzz ■ ■ ■ p-times. Thus, 



nf = o, 
n^ = 4a, 
n^ = o, zVi, 
nS = o, / / j. 

n^- = c*p, i ^ j, 
n& = o, i^j^k, 

n£U = 4«, i^jVi (io) 

Now, as discussed in [39] using H^ my n zP = U^ ynz p for all orders leads to some 
inconsistency in recovering the macroscopic fluid equations. To circumvent this 
issue, we use a factorized form for (central moment) attractors proposed in [39J. 
Essentially, in addition to satisfying Galilean invariance, the Maxwellian (equi- 
librium) satisfies the factorization property, i.e. independence along Cartesian 
coordinate directions, which immediately applies to its central moments. In 
the factorized central moment formulation, this property is extended to model 
non-equilibrium process, i.e. relaxation towards equilibrium. In other words, 
the higher order central moment attractors are given as its factorization in 
terms of lower order central moments that are not yet in equilibrium [39] . To 
proceed further, let us define the following post-collision continuous central 
moment of order (to + n + p): 

oo oo oo 

n 



2yn z p 



J J J m-u x r^ y -u y r^ z -u z rd^ x d^ z . (n) 



-oo — oo — oo 



Then, we consider the factorized form for attractors as 
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ft a * 


=n«n j = 0, 


frat 
LL ijk 


=fi i fi i S* = o, 


n at . . 

wo 


n j j n j j , 


frat 

iijk 




frat 

iijjk 


= UuTijjTik = 0, 


frat 

iijjkk 





Now, however, to correctly recover the momentum flux and pressure tensor 
in the macroscopic fluid dynamical equations, the diagonal components of 
the second-order central moments should preserve those obtained from the 
Maxwellian. That is, we set it"/ = c 2 s p. Thus, the 27 independent components 
of the local factorized central moment attractors can be written as 



frat 


= ft ai 


= fl at 

V 


= nf = 


ft ai 

XX 


= ft a * = 

yy 


fjai = 

22 


C 2 sP, 


xy 


= ft at = 

xz 


n a * = 

yz 


0, 


ft at 

xyy 


= U at = 

xzz 


-- fl at 

xxy 


= n a * 

yzz 



n; 
n 
ft 



at 

xxyy 



at 



at 



n 



yyzz 
at 



n 



xxyz 



at 

xyyzz 



n 



at 

xxyyzz 



n 



ai 



n 



at 



n 



at 

XJ/2 



n-XX-n-22 5 



nyyT\. zz , 



n 
ft 



at 

xyy2 



at 

xxyzz 



n 



at 

xyzz 



n 



at 

xxyyz 



0. 



(13) 



In essence, for the D3Q27 lattice, the fourth-order and sixth-order moments 
are factorized in terms of longitudinal second-order moments. It may be noted 
that symmetries in the factorization of the Maxwellian have been exploited to 
construct other types of quasi-equilibrium attractors recently [30] . 



Similarly for the continuous source central moments due to force fields, one 
possible choice is obtained by choosing that based on the local Maxwellian, 

P c; 



i.e. Af 1 



/ , which, however, leads to aliasing effects for higher 



order moments [40J. To circumvent this issue, a simple choice involves de- 
aliasing higher order moments while preserving its necessary effect on the 
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first-order central moments [40] which is extended to 3D in this work. Thus, 
we specify the continuous source central moments as 



m 



x myn z p 



1, n = 





P 


= 


0,n = 


1 


P 


= 


0,n = 





P 


= 1 



(14) 



Otherwise. 



4 Central Moment Lattice-Boltzmann Equation with Forcing Terms 



Let us now write the central moment lattice Boltzmann equation (LBE) with 
forcing terms by first defining a discrete distribution function supported by the 
discrete particle velocity set e a as f = \ f a ) = (fo, /i, /2, • • • , f2eY, a collision 
operator as f2 c = = (ft^, fij, fl^, . . . , f2§ 6 ) as well as a source term as 

S = \S a ) = (S , Si, S 2 , ■ ■ ■ , 5*26)^ based on central moments. The LBE can 
then be obtained as a discrete version of the continuous Boltzmann equation 
by temporally integrating along particle characteristics as jlD] 

t+i 

U-t + i> a ,t + i) = f a (-t, t) + n^ )t) + / s a{ ^ a6t+d) de. (is) 



In Eq. (15), the collision operator can be written in terms of the unknown 



collision kernel g projected to the orthogonal matrix of the moment basis 
as [36] 



J£ = n£(f,g) = (£•§)«, 



(16) 



where g = \g a ) = (#o;#i> <?2> • • • ,g2&)\ which will be derived later. The macro- 
scopic conserved moments, i.e. the local density and momentum, are obtained 
from the distribution function as 

26 

P=T,f*=(f*\P)> ( 17 ) 
26 

pUi = J2 fa^ai = (fa\e ai ) ,i G x,y,z. (18) 
We consider a semi-implicit representation for the source term, i.e. the last 



term in the above equation, Eq. (15), to provide second-order accuracy [IDJ, i.e. 
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( x + e a e,t+o) 



dO 



made effectively explicit by using the transformation f c 
it to 



S atf,t) + S atf + ^ a ,t + i)\ ■ This equation is then 

f a — oS a to reduce 



f a (t + i>*,t + i) = f a (-t, t) + n° ^ + 



(19) 



The explicit expressions for the source term S a as well as the collision ker- 
nel g will be derived so as to rigorously enforce Galilean invariance through 
a matching principle and a binomial transformation. These are discussed in 
Sees. [6] and 171 respectively. 



5 Various Discrete Central Moments and Galilean Invariance Match- 
ing Principle 



To facilitate the determination of the structure of the collision operator kernel 
g and the source terms S a , we now define the following discrete central mo- 
ments of the distribution function, Maxwellian, and source term, respectively: 







- u x 




~ u y, 


{&az 


- U z 






K at — 

lx x m y n zP 






ifiay 


~ u y, 


ifiaz 


- u z 


) P \fa), 




' = 




- u x 


\&ay 


~ u y] 


ifiaz 


- u z 


r\s a ). 


(20) 



Furthermore, the following definitions involving discrete central moments based 
on post-collision (f a ) and transformed (f a ) distribution functions, and its 

combination f a , are useful for further simplifications: 

K x m y n zP = \[e ax — u x ) {e ay — Uy) [e az — u z ) p \f a ) , 

K x m y n zV = \{e ax — u x ) m (e ay — u y ) n (e az — u z ) p \f a ) . (21) 

Based on the definition of the transformed distribution function as given in 
the last section, it immediately follows that 

K x m y n z p — K, x m y n z p — — U x m y n z p. (22) 



In order to preserve the main physical characteristic, i.e. Galilean invariance 
at the discrete level, we now invoke the key matching principle, which is to 
set the discrete central moments of the attractors of the distribution function 
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and the source terms, defined above, equal to their corresponding continuous 
central moments, whose forms are known exactly from the ansatz derived in 
Sec. [3j In other words, 



K x myn z p — TL x myn z p, (23) 

0~ x m y n zP ^ x rn y n zP ' (24) 

This yields the following expressions for the discrete local central moment 
attractors 
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(25) 



In addition, the discrete source central moments satisfy the following 



So = 0,5^ = F x ,a y = F y ,a z 



(7 x m y n zP — j 777f j Tbj J) ^> 1 . 



(26) 



Thus, finally, in view of Eq. (22), the attractors in terms of the transformed 
central moments can be written as 



^at „ ^at 1 ^at 1 j-, ^at 1 

K Q =[},K X = —-P x ,K y = —-P y ,K z = —-l* z , 
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^at ^at ^at ~ 

xyyzz xxyzz xxyyz w ' 

^xxyyzz ~ ^xx^yy^zz- (27) 



6 Various Discrete Raw Moments and Source Terms in Particle 
Velocity Space 



In order to construct an executable central moment LBM, the above informa- 
tion based on the central moments need to be related to the raw moments, i.e. 
those in the usual lattice or rest frame of reference. This can be readily accom- 
plished by means of the binomial theorem applied to the orthogonal products 
of the discrete quantities supported by the particle velocity set [40J. In this 
regard, the following notations that specify various discrete raw moments will 
be useful: 



f m n p l m n p \f\ foe 

x m y n zP J a^ax^ay^az \ c ax c ay c az\ J a I i \^ 



K 



V f e m e n e p = (e m e n e p I f ) (2V) 

/ j J a^ax^ay^az \^ax^ay^az\J a I ' \*" J J 



x m y n zP / j J a^ax^ay^az X^ax^ay^az 

! 

x m y n zP 



<T„ m „,n„v — S Q e. riT e. mi e v riz — (e Q , x e Q j / e^ z |S' Q ) . (30) 



Here and in what follows, the superscript "prime" (') is used to distinguish 
the raw moments from the central moments that are designated without the 



primes. Furthermore, similar to Eq. (22), the relation K x m v n = k 



x 'x m y n ~ '^x m y n zP 



\cr x m y n z p is satisfied. Based on the above, first, we write the raw moments 
of the distribution function of different orders supported by the particle veloc- 

ity set Tt x m y n zP = (f a \eax e ay e az) i n terms of the known quantities. To obtain a 
compact description of results, the following operator notation is helpful [3D] : 



a(f ai + fa 3 + fas + ■ ■ + Kf h +ffh+ffk + ---) + '-- 

= ( a E+ & E---W- (si) 

\ a a / 

where A = {an, a 2 , az 3 , • ■ •}, B = /3 2 , /3 3 , ■■■},■■ •. First, the conserved 
transformed raw moments follows directly from the definition as 



K = \fa\P) = P> K x= (faKx) = PU X ~ ^F x , 

«» = (f a \e ay ) = pu y - -F y , k z = (f a \e az ) = pu z - -F x . (32) 
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The non-conserved transformed raw moments can be written, using the above 



operator notation (Eq. (31)), in terms of the subsets of the particle velocity 



directions, which are presented in Appendix [B| 



The next step is to transform the central moments of the source terms (Eq. (26 )) 



in terms of raw moments by using the definitions, i.e. Eq. (20) and (30), which 
by the binomial theorem, readily yields the expressions that are enumerated in 
Appendix [Cj These moments should be related to the discrete source terms in 
particle velocity space so that an operational Galilean invariant approach can 
be derived. To accomplish this, we first obtain a set of intermediate quantities 
m%, which are the projections of the source terms to the orthogonal matrix 
of the moment basis /C, i.e. mp = (Kp\S a ), (3 = 0, 1, 2, . . . , 26, which can 
be obtained from the above using Eqs. (jlj) and (C.l). The details of are 
provided in Appendix |D} 

It is noted that m s g is equivalent to the following matrix formulation 



/C'S = (/C • S) a = ((K \S a ) , (KilSc) , (K 2 \S a ) , . . . , (K 26 \S a )) 

= (m s Q ,mlm s 2 , . . . ,m s 26 y = \m s a ) , (33) 

which can be exactly inverted by using the following orthogonal property 
of /C, i.e. /C _1 = K) ■ P -1 , where D is the diagonal matrix given by D = 
diag((K \K ) , (Ki\Ki) , (K 2 \K 2 ) , . . . , (K 26 \K 26 )) |H]. Exploiting this fact, the 
linear system (Eq. (33)) can be solved exactly to yield the expressions for 
the Galilean invariant source terms in velocity space S a in terms of m%, or 

equivalently the force ~F and velocity fields it. The final results of S a , where 
a = 0,l,2,...,26 are summarized in Appendix [El 



Finally, to obtain the collision kernel gp in the next section, we need to evaluate 
the expressions for its raw moments of various orders projected to the moment 
basis matrix JC, i.e. 

£(£ " iUZe n ay e p az = E (^ICCO h- ( 34 ) 

a p 



For conserved moments, it follows by definition that go = gi = g 2 = gz = 0. 
Again, exploiting the orthogonal property of /C, the moments of the collision 
kernel can be obtained which are presented in Appendix [F| 



The central moment LBE given in Eq. ( 19 ) can be rewritten in terms of the 



collision and streaming steps, respectively, as 

7 a &,t)=J a (?,t) + £T ^ + (35) 
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f a {-t + -? a ,t + l)=f a (t,t), (36) 

where the symbol "tilde" (~) in the first equation refers to the post-collision 
state. Furthermore, the conserved local fluid density and momentum are finally 
written in terms of the moments of the transformed distribution functions as 



26 _ _ 

p=E7 Q = (7», (37) 
26 - 1 - 1 

[>Ui=^2 f a e ai + -Fi = {f a \e m ) + -Fi, iex,y,z. (38) 

Q=0 — — 



7 Structure of the Central Moment Collision Operator 



We are now in a position to obtain the expressions for the collision kernel 
of the 3D central moment LBM in the presence of source terms. In essence, 
the procedure begins with the lowest order (i.e. second-order, off-diagonal) 

post-collision central moments (i.e. K xy ,K xz and K yz ), which are successively 



set equal to the corresponding attractors given in Eq. (|27|) (i.e. , k x \ and 



~Ky Z i respectively). This intermediate step is based on an equilibrium assump- 
tion. Dropping this modeling assumption to represent collision as a relaxation 
process by multiplying with a corresponding relaxation parameter results in 
the collision kernels g a for a given order [32] • Here, the relaxation parameter 
needs to be carefully applied to only those terms that are not yet in post- 
collision states, i.e. those that do not contain g@, where = 0, 1, 2, . . . , a — 1 
for a given g a . Then the results are transformed in terms of raw moments via 
the binomial theorem to yield expressions useful for computations. The details 
of various elements in obtaining the collision kernel are presented in jJU]. To 
simplify exposition, let us introduce the following notation: 

Vx m y n zP = K x m y n zP CT x m y n zPi (39) 



where the expressions for k 

X rn y n z'P 

and a„. are known from Sec. 6 In the 
following, for brevity, we present only the final results. For the above three 
off-diagonal central moments, we get 



9a 



95 



12 
12 



-Vxy + P u xUy + ^\p x Uy + a y u x ) 



1 

2' 



-r) xz + pu x u z + -{B x u z + a z u x ) 



(40) 
(41) 
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96 



12 



lyz 



pu y u z + - {d y ii z + a z u y 



(42) 



where 004, U5 and 00$ are relaxation parameters. Similarly, applying the pro- 
cedure for the remaining three second-order diagonal components with k xx = 
c 2 p, which preserves the Maxwellian values to provide the correct 



K yy — K zz 



momentum flux and pressure tensor, yields 



97 



UJ-j 

12 



~{Vxx - Vyy) + PK - V>1) + (K U X - PyUy) 



(43) 



98 = 



UJ8 

36 



(Vxx + Vyy ~ 2Vzz) + P( U l + «J ~ ^1) 



+(a x u x + a' y Uy - 2^m 2 ) 



(44) 



09 = 



18 



(Vxx + Vyy + Vzz) + p( U l + U l + U l) 

' z u z )+p . 



y^y 



(45) 



Next, carrying out the above matching, transformation, and relaxation steps 
(with the last of this applicable only for the pre-collision terms) successively 
to all the seven components of the third-order moments we get 



010 = 



^10 
24 



(Vxyy + Vxzz) + ^yVxy + UzVx*) + U AVyy + Vzz) 



-2pu x (u. 



~ 7}K( u l + u l) - u x {d' y u y + a z u z ) 



+ (u y g 4 + u z g 5 ) + -u x (-g 7 -g 8 + 2g 9 ), 



(46) 



0ii 



24 



-{v. 



xxy 



lyzzJ 



+ 2(UxVxy + U zVyz) + UyiVxx + Vzz) 



-2pu y (u 2 x + u 2 z ) - -v y (u 2 x + u 2 z ) - u y {a x u x + a z u z ) 
+ (u x g A + u z g 6 ) + ^u y (g 7 - g 8 + 2g 9 ), 



(47) 



912 



U12 

24 



(Vxxz + Vyyz) + ^xVxz + ^yVyz) + MVxx + Vyy) 
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~2pU z (u 2 x + lij) - 2 & *( U x + U l) ~ U z(K U X + OyUy) 

+ (Ux95 + Uyde) + T^Uzids + 9d) , 



(4£ 



#13 



W13 



(?7. 



■ry.i/ 



-2pu x {u 



xy 

1 / 2 2^ 

2^K 



.y;y 



u x {a y u y - cr z u z 



+3(u y g 4 - u z g 5 ) + -u x (-g 7 + 



(49) 



#14 



Wl4 



y 'lyzzJ 1 "\"-x'i X y "-z'lyzj ' ""y\ 

/ 2 
[11 

3 



-2pu y (ul - u 2 z ) - -a Jul - u*) - ^(5^ - a 2 < 



+3(«a;P4 - u z g 6 ) + -u y (g 7 + 3<? 8 ), 



(50) 



915 



^15 



- Vyyz) + 2 ( U xVxz ~ U yVyz) + U z(Vxx ~ Vyy) 



-2pu z (u 2 x - u 



/ 2 

3 



u z (cr x u x - a u 



y">VJ 



+3(u x g 5 - Uy%) + -u z g 7 , 



(51) 



916 



-»Ly* + M *?L + M 2/#^ + ^zTlxy ~ 2pU x U y U 



-^P'x u v u z + <?' y UxU z + a z u x Uy) 
3, ^ 

+ 2l M ^4 + Myfl'5 +U z g 6 ), 



(52) 



Notice that the cascaded structure is evident for the collision kernel starting 
from the third-order moments. Now, the next three diagonal components of 
the fourth-order central moments needs to carefully consider the non-zero 
factorized attractors given in terms of second-order components, i.e. K xxyy = 



K xx K yy: K xxzz ~ K xx K zzj anc i ^ yyzz 

corresponding collision kernels as 



KyyKzz ( see Eq. (27)). This yields the 
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_ Wi 7 

^ 17 " lT 



^/ ^/ ^/ / ^/ ^/ ^/ ^/ 
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(53) 
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(54) 
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xxyy 



Vxi 



+ 2 U X T] - U x 7] xzz + UyT) 



xxy 



-(ulVyy ~ ulVzz + «J^xx ~ U lVxx) ~ A ( U xUyVxy ~ U x U z Vxz) 
~\~{l^xx^yy K XX K ZZ ) + 3p(u x 1ly U X U Z ) 
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"V^y 
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6u x Uyg4 + Qu x u z g 5 



4 



+ -(-3u - 8)g 9 + 3u y gn - 3u z gi 2 + 2u x gi 3 + u y g 14 - u z g 15 . (55) 
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For calculating gn through gig in the above equations, we need the post col- 
lision states k xx , Kyy and k zz . These can be obtained from Eq. (22) as follows. 



^XX ^XX c^XXl 



K 



Ml ~ '' //'/ + .1° Mr 
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where the second-order transformed central moments, in turn, can be related 
to corresponding raw moments, which are known, as 



K>XX 



pu x - F x u x , 



— — — - 2 _ p 

— — — 2 TP 

K-zz — K zz pU z r Z U Z . 



Note that in terms of r] xm nzP these can also be written as 



K 3 



K 



yy ' 



Vxx + §97 + 6^8 + Q99 

Vyy ~ 697 + 6^8 + Q99 



- pu x - F X U X , 
~ P u y ~ ^y u yi 



Proceeding further for the remaining three fourth-order central moments using 

K xxyz = K xyyz = K xyzz = 0> We S e ^ 



920 = 



^20 



Vxxyz + U zVxxy + U yVxxz + ^ U xVxyz U y U zVxx ^ u xU z T] 



Ixy 



-2u x u y rj xz - u 2 x r) yz + 3pu\u y u z + a x u x u y u z + -u 2 x (a y u z + d' z u y ) 
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~o /3 , \^ 3 ^3 

-3u x u z g A - 3u x u y g 5 - I -u x + 1 1 g 6 - -u y u z g 7 - -u y u : 

3 3 3 1 1 

~^u y u z g 9 + -u z g u + -u y g 12 + -u z g 14 + -u y g 15 + 2u x g w , (56) 
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u x u z g s 



~ /3 , \ ^ ^3 ^3 

-3u y u z g A - \-u y + l\g 5 - 3u x u y g 6 + -u x u z g 7 - - 

3 3 3 1 1 

~^u x u z g 9 + -u z g 10 + -u x g 12 + -u z g 13 - -u x g 15 + 2u y g w , (57) 



21 



922 



^22 



Vxyzz + 2u zVxyz + U y f] xzz + U x f] U z f] 2u y U z T] i 
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3 ~ 3 ~ 3 ^ 1 ^ 1 ^ n ^ 
■-u x u y g 9 + -u y g w + -u x g n - -u y g 13 - -u x g u + 2u z g 16 . (58) 



The collision kernels for the three fifth-order central moments follow similarly 
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(61) 



Finally, for the one sixth-order component, we obtain the collision kernel based 



on the non-zero factorized attractor (see Eq. (27)) as 
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Note that the transformed raw moments of various orders, i.e. K x m y n zP and 
raw source moments, i.e. a xm „ zP needed for r] xm nzP for various m, n and 



p combinations can be obtained from Eqs. (32) and (B.l) and Eq. (C.l), 
respectively, which are given in Sec. |6j Similar to the 2D central moment 
LBM with source terms |4"0] . we can apply the Chapman- Enskog expansion 
to the above 3D formulation to show that its emergent dynamics corresponds 
to the Navier-Stokes equations representing fluid motion in the presence of 
general force fields. Some of the relaxation parameters in the collision model 
can be related to the transport coefficients. For example, those corresponding 
to the second-order moments control shear viscosity v of the fluid. That is, v = 
c « (w 57 — l) w here u u = ujj where j = 4, 5, 6, 7, 8. The rest of the parameters 
can be set either to 1 (i.e. equilibration) or adjusted independently to carefully 
control and improve numerical stability by means of a linear stability analysis, 
while all satisfying the usual bounds < cop < 2. 



8 Operational Steps of the Central Moment LBM 



To provide explicit expressions for the collision step in the central moment 



LBM as a stream-and-collide procedure (i.e. Eq. (35) and (36)), we first ex- 



pand the elements of the matrix multiplication of /C with g in Eq. (16). This 



yields the post-collision values of all the 27 components of the transformed 
distribution function in terms of the Galilean invariant collision kernel gp (see 



Sec. [7]) and source terms Sp (see Eq. (E.l) in Appendix. [E]) which can be 
summarized as follows: 



7o = 7o + [do - 2 99 + 4(?i 7 - 8(?26] + S , 
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(63) 



The above post-collision state allows completion of the streaming step via 



Eq. (36), following which frame-independent fields of 3D fluid motion can be 



obtained from Eqs. (37) and (38). In the implementation, various optimization 



strategies such as those discussed in [22J should be fully exploited to minimize 
the floating point operation count. 

Following the general outline of the above derivation, the central moment 
LBM was also formulated for the three-dimensional, fifteen velocity (D3Q15) 
lattice, which has a much reduced computational complexity when compared 



with the D3Q27 lattice. The results are summarized in Appendix G 



9 Numerical Tests 



Both the central moment formulations including forcing terms derived ear- 
lier, i.e. for the D3Q15 and D3Q27 lattices, were implemented and assessed. 
Let us now discuss the validation studies carried out for these computational 
approaches for a set of canonical problems for which analytical solutions are 
available. First, we consider a fully developed flow between parallel plates 
driven by a constant body force. The grid resolution was chosen to be 3 x 3 x 45 
with relaxation parameter u u = 1.818 for the second-order moments (u" = Uj 
where j = 4, 5, 6, 7, 8) that controls the kinematic viscosity u (= 0.0167 here). 
Rest of the relaxation parameters were set to be unity for this well 
as for all the simple canonical problems considered in our present numerical 
accuracy study. It may be noted that other values could be used for kinetic 
modes involving more complex situations (e.g. turbulent flows) and could also 
be optimized to improve numerical stability. For these parameters, three dif- 
ferent values of the body force were considered, i.e. F x = 2 x 10~ 7 ,4 x 10~ 7 
and 6 x 10 -7 corresponding to Reynolds numbers (based on the centerline ve- 
locity and half-width between plates) 3.6, 7.2 and 10.7, respectively. Half-way 
bounce back boundary condition was implemented to impose the no-slip con- 
dition at both the walls. Figure [2] shows a comparison between the computed 
results obtained using the central moment LBM implemented for D3Q15 and 
D3Q27 lattices with the analytical solution (u(z) = uq(1 — (z/L) 2 ), where L 
is the half-width and Uq = F x L 2 /(2v)). Excellent agreement is seen for both 
formulations with the benchmark analytical solution. Since the results with 
D3Q15 and D3Q27 lattices are essentially identical with the former involving 
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Fig. 2. Flow between parallel plates with a constant body force: Comparison of 
velocity profiles computed by D3Q15 (open symbols) and D3Q27 (filled circles) 
formulations of the central moment LBM with forcing term with analytical solution 
(lines) for different values of the body force F x . 

considerably lower operation count, henceforth we discuss the numerical per- 
formance only with the D3Q15 lattice. It may be noted that the advantage 
of the central moment formulation for this lattice, over the SRT approach 
lies in its enhanced numerical stability by independently and carefully adjust- 
ing the relaxation parameters for the kinetic modes. This and other assets 
such as better representation of kinetic layers are similar to the standard (raw 
moment) MRT approach. Comparison of such different collision models are 
subjects for future investigations. The central moment LBM using the D3Q15 
lattice was further assessed for the channel flow problem at higher Reynolds 
numbers. By considering the same resolution as before and setting the body 
force as F x = 1 x 10 -6 , two different Reynolds numbers of 111.8 and 447.2 
were considered by using u u = 1.923 and 1.961, respectively. Comparisons 
of computed and analytical solutions were made in Fig. |3j which again show 
good agreement. 

In order to quantify the error between the computed and analytical solutions 
and its variation at different resolutions, i.e. to establish the grid convergence 
of the 3D central moment LBM, the following test was carried out. We again 
considered channel flow with the computational domain discretized using 3 x 
3 x N nodes, where N is the number of nodes in the wall normal direction which 
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Fig. 3. Flow between parallel plates with a constant body force: Comparison of 
velocity profiles computed by D3Q15 formulation of the central moment LBM with 
forcing term (open symbols) with analytical solution (lines) for different values of 
Reynolds number Re. 

was varied. The parameters were chosen so as to satisfy diffusive scaling: the 
fluid velocity (or the Mach number) was made to scale with the resolution, i.e. 
Mo ~ Ax ~ 1/N. This ensures that the errors associated with compressibility 
effects also simultaneously reduce with increase in resolution. Thus, with a 
fixed viscosity v to maintain constant Reynolds number (Re = uqL/v) for 
different resolutions, using uq = F X L 2 / (2v) the body force scales as F x ~ 
1/L 3 ~ 1/iV 3 . Setting u u = 1.818 and considering F x = 6.958 x KT 6 for 
the coarsest resolution (N = 13) so that Re = 20.8, the relative errors in 
velocity field at different resolutions were computed using | \8u\ I2 = J2i II ( u c,i — 
U a,i)\\2/J2i \ \u a ,i\\2, where u c ^ and u a ,i are computed and analytical solutions, 
respectively, 1 1 - 1 1 2 is the standard second-norm and the subscript i represents 
the discrete location of the nodes. Figure [4] shows a logdog plot of the relative 
error as a function of the number of grid nodes. It is evident that quadratic 
grid convergence is maintained by the 3D cascaded LBM. 

We will now consider a different canonical problem, where the imposed body 
force is time dependent and thus represents a more stringent test of the central 
moment formulation derived in this work. In particular, flow between paral- 
lel plates driven by a force which varies sinusoidally in time was computed 
using this approach. If Q = 2it/T is the angular frequency, where T is the 
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Fig. 4. Grid convergence study of the D3Q15 formulation of the central moment 
LBM with forcing term for channel flow under diffusive scaling. Symbols represent 
relative (root-mean-square) error between the computed and analytical solution. 
Best fit slope of computed results is —1.96. 

time period of the application of the body force, it may be represented as 
F x = F m cos(Qt), where F m is its maximum amplitude. This problem, gen- 
erally termed as Womersley flow, is characterized by the dimensionless pa- 
rameter Wo = J~^L, also called as the Womersley number representing the 
relative effect of the unsteady response of the fluid flow to the imposed un- 
steady body force. It has the following analytical solution for the velocity 



profile u x (z) = 1Z 



iFm J 1 COS KH) 1 -tfit 



where (3 = y — iWo . Considering 



3 x 3 x 45 nodes and setting the maximum force amplitude F m = 1 x 10~ 5 with 
a time period of T = 10, 000, two different values of the relaxation parameter 
u u , i.e. 1.724 and 1.923, were used, which correspond to Wo of 3.3 and 6.6, 
respectively. Figures [5(a) | and 5(b) show comparisons of the computed veloc- 



ity profiles with the above analytical solution for different instants within the 
first half of the time period T at these two Womersley numbers. It is clear 
that the central moment LBM reproduces the sharp variations in the velocity 
profiles at different instants as prescribed by the analytical solution, with very 
good agreement found between them. Furthermore, the variations in both the 
amplitude as well as the lag of the response of the fluid flow as seen by its 
velocity profiles at different Womersley numbers are well reproduced by the 
computational approach presented in this work. 
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It may be noted that in all the problems considered above, the velocity field 
has variation along only one coordinate direction normal to the direction of 
the driving body force. Thus, as a final example, we consider fully developed 
flow through a square duct in which the flow field has variations in both the 
coordinate directions normal to the direction of application of the driving 
force. It has the following analytical solution for the velocity field given in 
terms of an infinite orthogonal (Fourier) series 



cosh ( 



(2n-l)nz 
2a 



cosh ( 



(2n-l)7T 
2 



COS 



^ (2n-l)ny ^ 

(2n-l) 3 



(64) 



where —a < y < a and —a < z < a. Here, a is the duct half-width. We 
considered the square duct to be resolved by using 3 x 45 x 45 nodes. A 
constant body force of magnitude F x = 1 x 10~ 6 was applied by considering 
the relaxation parameter uf equal to 1.923 such that the Reynolds number 
(based on maximum or centerline velocity and duct half-width) is equal to 
65.7. As before, the no-slip condition at the walls was imposed using the half- 
way bounce back approach. Figures 6(a) and |6(b) show a comparison between 
the surface contours of the computed and analytical solution of the velocity 
field for the above condition. It is seen that the 3D central moment LBM with 
forcing term is able to reproduce the distribution of the velocity field over the 
cross-section of the square duct. In order to more clearly make a quantitative 
comparison, Fig. [7] shows plots of the computed velocity profiles at different 
locations in the cross-section of the duct and their comparison with the cor- 
responding analytical solution (see Eq. (64)) Evidently, the results computed 
using the central moment LBM are found to be in excellent agreement with 
the benchmark solution. 



10 Summary and Conclusions 



A derivation of the 3D central moment lattice Boltzmann method (LBM) in 
the presence of forcing terms is presented. Suitable orthogonal moment ba- 
sis for the D3Q27 and D3Q15 lattices are chosen for the specification of the 
local attractors and source terms in terms of central moments. In particular, 
recently proposed factorized form of local attractors for higher moments and 
de-aliased source terms that influence only conserved moments, which are ob- 
tained from modifications of the properties of the Maxwellian are considered in 
the construction of the approach. A Galilean invariance matching principle is 
invoked that exactly preserves the continuous central moments of the attractor 
and the source terms at the discrete level for all orders supported by the par- 
ticle velocity set. Based on these, expressions for the temporally semi-implicit 
and second-order accurate sources are derived through an exact inversion due 
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to the orthogonal properties of the moment basis. The central moment LBM, 
whose elements are equivalently expressed in terms of raw moments using the 
binominal theorem, represents frame independent fluid motion in the presence 
of general external or self-consistent internal forces. A set of numerical tests 
was carried out for problems involving channel flow driven by constant and 
temporally varying (periodic) body forces, and flow through a square duct to 
assess the accuracy of the central moment LBM with forcing term derived in 
this paper. It is demonstrated that the method maintains second-order grid 
convergence under diffusive scaling. Comparisons of the computed results are 
found to be in excellent agreement with analytical solutions for all the bench- 
mark problems considered. 
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(b) 

Fig. 5. Flow between parallel plates with a temporally varying body force: Compar- 
ison of velocity profiles computed by the D3Q15 formulation of the central moment 
LBM with forcing term (open symbols) with analytical solution (lines) at different 
instants within a time period T. (a) Wo = 3.3 and (b) Wo = 6.6, where Wo is the 
Womersley number. 32 




(a) 



Analytical 




(b) 



Fig. 6. Flow through a square duct with side length 2a subjected to a constant body 
force: Comparison of surface contours of the velocity field for Reynolds number 
Re = 65.7 (a) computed by the D3Q15 formulation of the central moment LBM 



with forcing term with (b) analytical solution (see Eq. (64)) 
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Fig. 7. Flow through a square duct with side length 2a subjected to a constant 
body force: Comparison of velocity profiles computed by the D3Q15 formulation 
of the central moment LBM with forcing term (symbols) with analytical solution 
(lines) (see Eq. (64)) at different locations in the duct cross-section for Reynolds 
number Re = 65.7. 
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A Appendix: Orthogonal Matrix of the Moment Basis /C for the 
D3Q27 Lattice 

A main element of the central moment method is the moment basis. The 
components of the orthogonal matrix of the the moment basis derived in Sec. [2] 
(see Eq. Q) can be written as 
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B Appendix: Non-conserved Transformed Raw Moments for the 
D3Q27 Lattice 

The non-conserved transformed raw moments of various orders are given in 
terms of the subsets of the particle velocity directions as 
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E7^ 

Q = 

26 



a a 



e ay e az 



a a 
'A 14 B 14 , 



a a 



(f a\ e ay e az) E fa e ay e az 



o=0 



f ai 



A 5 S 5 

E-E]®7 a , 



E-E]®7o, 



E-E)®7o, 

a a / 

An Bli\ 

E-E ®7a, 



(/al e aa; e o?/) E f a e ax e ay ( E E I ® /a> 



^13 Bi 3 '' 

E-E] ®/«> 



E-E ®/a, 



A15 B 1S > 
E-E] ®/«> 



26 



(/al e aa; e a?/Ca2) /a 



Ai 6 Bi 6 N 

= (E-E]®7 



a=0 



a a 



26 



a=0 
26 

E7 

Q = 

26 



2 „2 
ay 



11 | e 2 e M = Vf e 2 e 

(,/al e a:r e az) 



2 2 
P P 
a otx az 



EJ ® /«> 

Ai 8 \ 

E ®/«. 



(/al e aj/ e az/ 



'Aig 



E7,ei.e? 



o=0 



a ay az 



E ®/«. 
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- 2 ^ - 2 

K xxyz = (f a\ e ax e ay e az) = ^ f a ^ ax e ay e az 

Q = 



26 

(/al e «a; e aj/ e a2) = / a e »:c e m/ e a2 



a=0 

-' - 2 ^ " 2 

a=0 



2 2 26 - 2 2 

a=0 




K xxyzz = (f a\ e ax e ay e az) = f a e ax e ay e az = I ~~ I ® /a' 

a=0 \ a a / 

_ 26 _ /A 25 B 25 \ _ 

(f a\ e ax e ay e az) = f a e ax e ay e az = I _ 51 I ® ^a' 

a=0 \ a a J 

^xxyyzz ~ (fa\ e ax e ay e az) ~ f a e ax e ay e az ~ I _ I ® /a' (^-1) 

a=0 V a a I 



xxyyz 



where 



A 4 = {7, 10, 19, 22, 23, 26} , B 4 = {8, 9, 20, 21, 24, 25} , 
A 5 = {11, 14, 19, 21, 24, 26} , B 5 = {12, 13, 20, 22, 23, 25} , 
A 6 = {15, 18, 19, 20, 25, 26} , B 6 = {16, 17, 21, 22, 23, 24} , 
A 7 = {1, 2, 7, 8, 9, 10, 11, 12, 13, 14, 19, 20, 21, 22, 23, 24, 25, 26} , 
A 8 = {3, 4, 7, 8, 9, 10, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26} , 
A 9 = {5, 6, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26} , 
A w = {7, 9, 19, 21, 23, 25} , B w = {8, 10, 20, 22, 24, 26} , 
A u = {11, 13, 19, 21, 23, 25} , B u = {12, 14, 20, 22, 24, 26} , 
A 12 = {7, 8, 19, 20, 23, 24} , B 12 = {9, 10, 21, 22, 25, 26} , 
A 13 = {15, 17, 19, 20, 23, 24} , B 13 = {16, 18, 21, 22, 25, 26} , 
A u = {11, 12, 19, 20, 21, 22} , B u = {13, 14, 23, 24, 25, 26} , 
A 15 = {15, 16, 19, 20, 21, 22} , B 15 = {17, 18, 23, 24, 25, 26} , 
A 16 = {19, 22, 24, 25} , B w = {20, 21, 23, 26} , 
A 17 = {7, 8, 9, 10, 19, 20, 21, 22, 23, 24, 25, 26} , 
A 18 = {11, 12, 13, 14, 19, 20, 21, 22, 23, 24, 25, 26} , 
A 19 = {15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26} , 
A 20 = {19, 20, 25, 26} , B 20 = {21, 22, 23, 24} , 
A 21 = {19, 21, 24, 26} , B 21 = {20, 22, 23, 25} , 
A 22 = {19, 22, 23, 26} , B 22 = {20, 21, 24, 25} , 
A 23 = {19, 21, 23, 25} , B 23 = {20, 22, 24, 26} , 
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A 24 = 
A 25 = 
A 26 = 



{19, 20, 23, 24} , B 2A = {21, 22, 25, 26} , 
{19, 20, 21, 22} , B 25 = {23, 24, 25, 26} , 
{19,20,21,22,23,24,25,26}. 



C Appendix: Raw Source Moments for the D3Q27 Lattice 



The raw source moments of various orders are given in terms of the Cartesian 
components of the force field as 



v'o = (S a \p) = 0, 

®x (Sa\&ax) F x , 

Uy = (S a \e ay ) = Fy, 

O 'z i^a\&az) F Z) 

a xx = (S a \e 2 ax ) = 2F x u x , 
o,yy = (S a \e ay ) = 2F y u y , 
Kz = (S a \e 2 az ) = 2F z u z , 

® xy ($ot \&ax&ay) F x lly -\- FyU x , 
® xz {Sa\&ax&az) F X U Z -\- F Z U X) 
Cy Z (S a | G-ay&az) FyU z -\- F z Uy, 

^xyy = (Sa\ e ax e ay) = F x U y + 2F y U y U x , 

Kzz = (S a \e ax e 2 az ) = F x u\ + 2F z u z u x , 

®xxy ~ Wa\ e ax e oty) — F y U x 2F x U x U y , 
®yzz = (^a\ e ay e az) = FyU z + 2F z U z U y , 
a xxz = (Sa\ e ax e ay) = F Z U X + 2F X U X U Z , 
® yyz (^q I ^ay^otz) F z Uy -\- 2FyUyU z , 

~ xyz (S a \6. ax G a yG az ) F x UyU z -\- FyU x tl z -\- F z U X Uy, 
®xxyy = (Sa\ e ax e ay) = ^F x U x U y + 2F y U y U x , 
a xxzz = (Sa\&ax e az) = ^F X U X U Z + 2F Z U Z U X , 
a yyzz = {^o\ e ay e <xz) = ^F y U y U z + 2F z U z U y , 
a xxyz = (Sa\e a x e ay e az) — u x (F y U z + F z U y ) + 2F x U x U y U z , 
® xyyz (S a \& ax 6 a y6 az ) Uy(F x 11 z -\- F z Ux) ~\~ 2FyUyU x U z , 
' xyzz {^a\^ax^ay^ az ) ^zi-Fx^y FyU x ) 4~ 2F z U z U x lly, 
®xyyzz = (^a\ e ax e a y e az) = F x U y U z + 2u x UyU z (F y U z + F z Uy) , 

Kxyzz = ( s a\e 2 ax e ay e 2 az ) = F y u 2 x u 2 z + 2u x u y u z (F x u z + F z u x ), 
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'xxyyz ~ (Sa\& ax e ay e az) — F z U x U y + 2u x U y U z (F x U y + F y U x ), 

®xxyyzz (Sal^ax^ay^az) 2u x UyU z {F x UyU z + FyU z U x 4~ F z U x Uy) . (C.l) 



D Appendix: Projections of the Raw Source Moments to the Or- 
thogonal Basis Vectors for the D3Q27 Lattice 



The orthogonal projections of the raw source moments can be written as 



m s = 


~- (K 


S a 


m\ = 


-- (K 




m s 2 = 


~-{K 2 


S a 


ml = 


--(Ks 


S a 


m\ = 




S a 


ml = 


~-(K 5 


S a 


ml = 


--{K 6 




m s 7 = 


--(K 7 


S a 


ml = 


~-(Ks 


S a 


ml = 


--(K 9 


S a 


m s 10 = 


1 r/ 

(#io 


S a 


" b ii 


(Km 
\-* v 11 


s„ 

'-'a 


m 12 = 


(K\2 


Q 

o a 


m s 13 = 


(K 13 


S a 


m s 14 = 


(K u 


S a 




(K 15 


S a 


m s 16 = 


(K m 


S a 


m{ 7 = 


(K 17 


S a 


m s 18 = 


(K 18 


S a 


m s 19 = 


(K w 


S a 


™20 = 


(K 20 


S a 


™>21 = 


(K 21 


S a 


™ s 22 = 


(K 22 


S a 



yi 



(F x Uy ~\~ FyU x ), 

= (F x u z + F z u x ), 

= (FyU Z + F Z Uy), 

^yFxIlx FyUy), 
= 2(F x u x + FyUy - 2F z u z ), 
= 2(F x u x + FyUy + F z u z ), 

= 3 [F x {ul + u 2 z ) + 2u x [F y Uy + F z u z ) 
= 3 \F y (u 2 x + u 2 ) + 2u y {F x u x + F z u z ) 
F z (u 2 x + u 2 y ) + 2u z (F x u x + FyUy) 
F x (u 2 y - u 2 z ) + 2u x (F y Uy - F z u z ) 
F y( u l ~ u l) + 1u y (F x u x - F z u z ) 
F z (u 2 x - v 2 ) + 2u z {F x u x - FyUy) 

-^x^y^z ^y^x^z ^ z^x^ 



-*F X , 
-4F Z , 



F x u x (u 2 y + u\) + F y Uy(u 2 x + u 2 z ) + F z u z (u 2 x + u 2 y ) 

8(F X U X + FyUy + F z u z ), 

F x u x {u 2 y + u\) + F y Uy(u 2 x - 2u 2 z ) + F z u z (u 2 x - 2u 2 y ) 
A(2F x u x FyUy F z u z ), 



= 6 



6 



F x u x (u 



+ U l(FyU y - F Z U Z ) 



F z u z ), 



-A{F y u y 

(3u 2 x - 2) [F y u z + F z u y ] + 6F x u x u y u z , 
(3u 2 y - 2) [F x u z + F z u x \ + 6FyU y u x u z , 
(3u 2 z - 2) [F x u y + F y u x ] + 6F z u z u x u y , 
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rn 



25 



m 



26 



(K 26 \S a ) 



m s 23 = (K 23 \S a ) = 9R 

+18u x 
m s 24 = (K 24 \S a } = 9F y l 



2 2 
U y U z 



y y + Ul) 



FyUyUl + F z U Z Uy - ^{FyUy + F, 



1 2 

u x u z 



u , 



+18% 



F x u x u\ + F z u z ul 



(K 25 \S a ) = 9F z 



2 2 ^ / / 2 , 2\ 

u x u y ~ 3 I K + V 



+18w 



2 

3 
2 

3 
2 

3 
2 

2 

3 
2 



F X U X Uy + FyUyU x (F X U X + FyUy^ 



5Au 2 y ui - 36(w; + <) + 24 



+F y % 54«X - 36 (< + wf) + 24 



54m>; - 36(< + iiy) + 24 



(D.l) 



E Appendix: Source Terms in Particle Velocity Space for the D3Q27 
Lattice 



The source terms in particle velocity space obtained by solving Eq. (33) are 
given by 



50 = 216 [8 ™° ~ ^ + 24 ™ Sl7 ~ ^ > 

51 = 216 + 12 ™ &1 + 18m&7 + 6 ™ 8 ~ 12m&9 ~ 12 ^ Sl ° ~ 12 ^ 8 

+ 12m* 3 + 4m* 6 ] , 

52 = 21~6 ^ 8 ™ S ° ~ + + 6 ™ 8 _ 12 ^ + 12 ^ Sl0 ~~ 12 ^ 8 

-12m^ 3 + 4my , 

^ 3 = 2~l6 ^ 8 ™ S ° + 12 ™ &2 ~ 18 ™ Sj + 6 ™ 8 ~ 12 ^ ~ 12 ^ x + 6 ^" 8 



-18m* 9 + 12m* 4 + 4m* 6 ] , 
1 
2T6 

^19 — J-^'"-24 T '±'/t 2 6j ) 



S A = — [8m* - 12m* - 18m* + 6m* - 12m*, + 12m^ + 6m* 18 
-18mi 9 - 12m*, 4 + 4m?, 

Sb = 216 + 12 ™ &3 ~ 12m&s ~ 12 ™ &9 ~ 12m&12 + 6 ^ 8 + 18 ^" 9 



+ 12m* 5 + 4m* 6 ] , 
1 

2T6 



S 6 = [8m* - 12m* - 12m* - 12m* + 12m* 12 + 6m* 8 + 18m* 19 
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10 



'21 



-12m s 25 + 4m s 26 \ , 

[8m s + 12m* + 12m* + 18m* + 12m* - 3m* - 3m* x 

216 

+27m* 3 + 27mi 4 - 6m* 7 + 3m 48 + 9m* 9 - 18m 22 - 6m 23 
-6m 24 - 2m 26 ] , 

[8m 9 - 12m* + 12m 2 - 18ml + 12m 8 + 3m* - Sm? n 

-27m* 3 + 27m* 4 - 6m* 7 + 3m* 8 + 9m* 9 + 18m* 2 + 6m 23 
-6m 24 - 2m 26 ] , 

-j- [8m* + 12m* - 12m* - 18m* + 12m* - 3m* + 3m* x 

+27m* 3 - 27m* 4 - 6m* 7 + 3m* 8 + 9m* 9 + 18m* 2 - 6m s 23 
+6m* 4 - 2m* 6 ] , 

^ [8m* - 12m* - 12m* + 18m* + 12m* + 3m* 10 + 3m* x 

-27m^ 3 - 27m 44 - 6m^ 7 + 3m* 8 + 9m^ 9 - 18m 22 + 6m 23 
+6m* 4 - 2m* 6 ] , 

-j- [8m 9 + 12m* + 12m? + 18m? + 18m 7 - 6m 8 - 3m 
216 

-3m* 2 - 27m^ 3 + 27mf 5 - Qm s 17 + 3m s ls - 9m s ig - 18m 
-6m 23 - 6m 25 - 2m 26 ] , 

[8m* - 12m* + 12m* - 18m* + 18m*, - 6m* + 3m* 

-3m^ 2 + 27mf 3 + 27mf 5 - 6m s 17 + 3mf 8 - 9m s 19 + 18m*. 
+6m 23 - 6m 25 - 2m 26 ] , 

^ [8m* + 12m* - 12m* - 18m* + 18m* - 6m* - 3m 

+3m 42 - 27m* 3 - 27m* 5 - 6m* 7 + 3mi 8 - 9m* 9 + 18m 
-6m* 3 + 6m* 5 - 2m* 6 ] , 

^ [8m* - 12m* - 12m* + 18m* + 18m* - 6m* + 3m 

+3m 42 + 27m 43 - 27m* 5 - 6m 47 + 3m* 8 - 9m* 9 - 18m 21 
+6m* 3 + 6m* 5 - 2m* 6 ] , 

[8m* + 12m* + 12m* + 18m* - 18m* - 6m* - 3m* x 

-3m* 2 - 27m* 4 - 27mi 5 - 6m* 7 - 6m* 8 - 18m 20 - 6m 24 
-6m* 5 - 2m* 6 ] , 

-J- [8mg - 12m*. + 12m*. - 18m*; - 18m*. - 6m s 8 + 3m^ 

-3m* 2 + 27m* 4 - 27m* 5 - 6m* 7 - 6m 48 + 18m 20 + 6m 24 
-6m s 25 - 2m s 26 ] , 

-j- [8m* + 12m* - 12m* - 18m* - 18m* - 6m* - 3m* x 
216 



21 



10 
21 



10 



41 



+3m s 12 - 27m* 4 + 27m* 5 - 6m s 17 - 6m s ls + 18m 20 - 6m 24 
+6m* 5 - 2m* 6 ] , 



[8m* ~ 12m* - 12m* + 18m* - 18m* - 6m* + 3m* x 

+3m s 12 + 27m s u + 27m* 5 - <6m\ 7 - 6m* 8 - 18m 20 + 6m 24 
+6m* 5 - 2m* 6 ] , 



[8m* + 12 (m* + m* + m*) + 18 (m* + m* + m*) 

+ 12mg + 6(m* + m*^ + m* 2 ) + 27m^ 6 + 6m* 7 + 9(m 20 
+m 21 + m 22 ) + 3(m 23 + m 24 + m 25 ) + m 26 ] , 

-J- [8m* + 12(-m* + m* + m*) + 18(-m* - m* + m*) 

+ 12m* + 6(-m* + m* x + m* 2 ) - 27m* 6 + 6m* 7 + 9(m*, 
-m s 21 - m s 22 ) + 3(-m*, 3 + m^ + m| 5 ) + m s 2e ] , 

1 [8m*, + 12(m* - m 2 + m 3 ) + 18(-m4 + mjj - m@) 



+ 12m*, + 6(m^ - m s n + m\ 2 ) - 27mf 6 + 6m s 17 + 9(-m 20 
+m 21 - m 22 ) + 3(m 23 - m 24 + m 25 ) + m 26 ] , 

-j- [8m* + 12(-m* - m* + m*) + 18(m* - m* - m*) 
21b 

+ 12mg + 6(-m* - m^ + m\ 2 ) + 27m* 6 + 6m s 17 + 9(-m 20 
-m 21 + my + 3(-m 23 - m 24 + m 25 ) + m 26 ] , 

-j- [8m* + 12(m* + m* - m*) + 18(m* - m* - m*) 
21b 

+12mg + 6(m* + m^ - m* 2 ) - 27mi 6 + 6m^ 7 + 9(-m 20 
-m^ + m* 2 ) + 3(m* 3 + m* 4 - m* 5 ) + m* 6 ] , 

-j- [8m* + 12(-m* + m* - m*) + 18(-m* + m* - m*) 

+ 12m* + 6(-m* + m* x - m* 2 ) + 27m* 6 + 6m* 7 + 9(-m* 
+m 21 - m 22 ) + 3(-m 23 + m 24 - m 25 ) + m 26 ] , 

-j- [8m* + 12(m* - m* - m*) + 18(-m* - m* + m*) 
21b 

+ 12m*, + 6(m* - m s u - m{ 2 ) + 27mi 6 + 6m* 7 + 9(m 20 
-m s 2l - m* 2 ) + 3(m* 3 - m* 4 - m* 5 ) + m* 6 ] , 

-j— [8mo + 12(-m* - m 2 - mg) + 18(m 4 + mg + m|) 

+12mg + 6(-m s w - m s u - m* 2 ) - 27m* 6 + 6m^ 7 + 9(m 20 
+m 21 + m 22 ) + 3(-m 23 - m 24 - m 25 ) + m 26 ] . 
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F Appendix: Moments of the Collision Kernel for the D3Q27 Lat- 
tice 



The moments of the collision kernel follow from the orthogonal property of 
the moment basis matrix /C, which are given by 



E(/c-g) Q = E(^lp)^ : 

a 3 

Y (Kp\e ax )g~p- 



0, 

o, 

5Z(/C • g) a e ay = Y ( K p\e*y) 9/3 = 0, 

0, 



Y^ ' §) Y ( K pKz)g/3 

a p 

£(/C ■ g) a e ax e ay = Y ( K i3\eaxe ay ) dp = 12<? 4 , 

a p 

£(/C • g) a e ax e az = Y (Kp\e ax e az ) gp = 12g 5 , 

a 3 

£(/C ■ g) a e ay e az = Y (Kp\e ay e az ) gp = I2g 6 , 

a p 

£(/C ■ g) a e 2 ax = Y ( K p\ e lx) gfi = 6^7 + 6^ 8 + 6g 9 , 

a p 

• s) a e 2 ay = Y ( K p\ e l y ) dp = - Q gr + 6^s + 6^9, 



a p 

£(K-g) Q eL = £(i^|eL>^ 

a p 



£(/C • g) a e ax e 2 ay = Y (Kp 

a p 

]T(/C • g) a e ax e 2 az = Y ( K P 

a p 

Y(K ■ i)y ax e ay = Y ( K P 

a p 

Yi K ' S)ae ay e 2 az = Y ( K P 

a p 

• g)*e 2 ax e az = Y ( K P 

a p 

£(£ ■ S) a e 2 ay e az = Y ( K P 



e ax e ay ) dp : 



)h = 



€-a X €- az 



e a X e ay) gp : 



e e 2 



elxCaz) dp : 

e 2 ay e az ) gp ■■ 



^ ' (^-~ ' g)a^a X &ay&az ^ ] p\& ax & ay & c 



]T(/C • g) a e 2 ax e 2 ay = Y ( K p\ e lx e l y ) dp 

a p 



-12<7 8 + 6^9, 
12^io + 4^i3, 
12^io - 4^i3, 
12g u + 4^4, 
12g n - 4g u , 
12g 12 + 4g 15 , 
12g 12 - 4<7i 5 , 

8<7l6, 

8^ 8 + 8^9 + 4^i 7 
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+4<7i8 + Ag X g, 

E(£ • g) a e 2 ax e 2 az = E (Ki3\e 2 ax e 2 az ) gp = 4g 7 - 4g s + 8g 9 

+4g 17 + Ag 18 - Ag 19 , 
E(^ ' S) a e 2 ay e 2 az = E ( K p\ e l y e lz) h = ~ 4 97 ~ 4^8 + 8# 9 



a 


P 

g;\ 2 \ / jr 


G/yil 7 / 0/-? 


+4^17 - 8^i 8 , 
8 Or + 8(720, 


E^ 

a 


^ax^ay^-az/ 9f3 


8fi- 5 + 85-21, 


£(£ 

a 


§) a^ax^ay^az ^ ] (-^/? 


^ax^ay^az) 9)3 


8^4 + 8^22, 


a 




e ax&ay e az) 9/3 = 


16^io + 8^23, 


£(£ 

a 




e ax e ay e az) 9/3 = 


16^11 + 8^24, 


£(£ 

a 


i) a e 2 ax e 2 ay e az = E 

g)aeL e aj/ e L = E 


e ax e ay e az) 9/3 = 


16^12 + 8^25, 


£(£ 


2 2 2 \ - 
^ax e ay e az) 9/3 


8^9 + 8^17 + 8£ 2 6. (F.l) 



G Appendix: Formulation of the Central Moment LBM for the 
Three-dimensional, Fifteen Velocity (D3Q15) Lattice 



G.l Moment Basis 



The particle velocity for the D3Q15 lattice ~~t a (see Fig. G.l) is given by 



wv 



(0,0,0), a = 

(±1,0,0),(0,±1,0),(0,0,±1), a = l,---,6 
(±1,±1,±1), 



a = 7, • • • , 14 



(G.l) 



The components of the nominal moment basis chosen are 



\To) = 


\p) = \[ta\ Q ) 


\Ti) = 




\Ti) = 




\T,) = 
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Fig. G.l. Three-dimensional, fifteen particle velocity (D3Q27) lattice. 



IT 4 ) 




^ctx^-ay) i 


\T 5 ) 




^ctx^-az) i 


\T 6 ) 






\Ti) 




2 2 \ 
P — P ) 

^ax ^ayl i 


\Ts) 




2 2 \ 
P — P ) 

^ax ^azl i 


\T 9 ) 




2 _|_ 2 _|_ 2 \ 

^ax ' ^ay ' ^azl ) 


\Tio) 




e ax( e ax e ay e az)) > 


\Tn) 




e ay(&ax e ay e az)) ? 


\Tn) 




e f e 2 +e 2 + e 2 ^ 


\Tis) 




dax&ay&az) j 


\Tu) 




2 2,2 2,2 

e ax e ay ' e ax e az ' e ay e 



Based on the above set, the components of the orthogonal basis vectors are 
obtained by means of the Gram-Schmidt procedure, which are given by 



K ) = \p) = \\-t a \°), 
Ki) = \e ax ) , 
K 2 ) = \e ay ) , 
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| ^OLZ j t 




\Ka) 






\Kf>) 






\K 6 ) 


l&ay&az) i 




\K 7 ) 


1 2 2 \ 




\K%) 


| 2 | 2 i 2 \ o|2\ 
\ e ax ' e ay ' ^azl " \ e az) i 




\Kv) 


= \ e ax + e ay + e az) ~ 2 \p) , 




\K W ) 


= 5|eax(eL + e^ + eL))-13 




\Kii) 


= 5 e Qy (e Q;r + e a?/ + e az j) — 13 




\Kl2) 


= 5 e a2 (e as + e ay + e az )) — 13 




\Kn) 


l&ax&ay&az) j 




\Ku) 


qj"i 1 2 2 1 2 2 1 2 2 \ 

— ou e Qa ,e ay + e aa; e Q2 -+- e a y e az) 


"40|eL 



,2 
'ay 



They can be written in a corresponding matrix form as /C 
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where 



K = [\K ) , \K,) , \K 2 ) , \K 3 ) , \K 4 ) , \K 5 ) , \K 6 ) , \K 7 ) , \K 8 ) 
\K 9 ) , \K 10 ) , \K n ) , \K 12 ) , \K 13 ) , \K U )\ . 



(G.3) 



G.2 Various Raw Moments and Source Terms in Velocity Space 



The above orthogonal matrix results in a set of moments of the collision kernel, 
which are needed in the construction of the collision operator, and are given 
by 
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E(£ • g) a e a:r = ^ (i^e^) = 0, 

a p 

E(^ ■ g)a e «y = E (^l e «s/) & = °» 
E(£ • g) Q e az = ^ (^|e az ) = 0, 

a p 

E(£ ' S)ae ax e ay = E (-^le^e^) gp = 8g 4 , 

a p 

E(£ ■ g) a e ax e az = (Kp\e ax e az ) gp = 8g 5 , 

a p 

E(£ • S)ae ay e az = (Kp\e ay e az ) gp = 8g 6 , 

a p 

£(/C • g) a e 2 ax = E (Kp\e 2 ax ) gp = 2g 7 + 2g 8 + Qg 9 , 

a p 

E(£ • S) a e 2 ay = E ( K p\ e ly) 9p = -297 + 2(fe + 6# 9 , 
£(/C • §)«£, = E (^|eL) & = -4^ 8 + 6g 9 , 



E(^ • S) a e ax e 2 az = ( K P 

a p 

E(^ • S) a e 2 ax e ay = ]T (Kp 

a p 

E(^ ' S) a e ay e 2 az = E ( K P 

a p 

E(^ • S) a e 2 ax e az = ( K P 

a p 

E(^ • S) a e 2 ay e az = £ ( K P 

a p 

^ j (J^ ' &)aC ax 6 ay 6 az ^ ] (Kp\Cc& 
a p 

E(^ ' g%eLe^ = E 

a p 

£(/C • g) a eLeL = E <*> 
E(* • S)«e 2 ay e 2 az = E 



e a xel y )gp = 16^io, 

e ay e L)^ = 16 ^ll> 
eLe a2 )^ = 16^12, 
e 2 ay e az )gp^l6g 12: 
e ay G az ) gp = 8gis, 
e 2 ax e 2 ay )gp = 8g 9 + 16<?i 4 , 
e 2 ax e 2 az )gp = 8g 9 + 16g u , 



e 2 ay e 2 az )gp^8g 9 + 16<?i4. 



(G.4) 



Note that unlike the D3Q27 lattice, additional degeneracies for various third 
and higher order moment basis vectors exist for the D3Q15 lattice, as it con- 
tains a more limited set of independent basis vectors. 
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It may be noted that the components of the raw moments of the source terms 
<7 x .m y n 2 p = (Sa\ e ax e ay e az) due ^° f° rce fields can be obtained in an analogous 
manner as determined for the D3Q27 lattice (see AppendixjC]). The projections 
of the source terms to the orthogonal matrix of the moment basis /C, i.e. 
(K p \S a ), = 0, 1, 2, . . . , 14 for this lattice yield 



m c 
ml 
m\ 
m\ 
m\ 
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m s - t 
m\ 
ml 
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13 



ni 



ii 



--(K 
--{Ki 
~-{K 2 
~-(K 3 
-- (K 4 
'-{K 5 
--(K 6 
--(K 7 
~-(K 8 
~-(K 9 

(K w 
(Ku 
(K 12 
(K 13 
(K 14 



S a 
S a 
S a 
S a 
S a 
S a 
S a 
S a 
S a 
S a 

S a 

S a 

S a 
S a 
S n 



y^x ) i 



(F x u z + F z u x ), 

(FyU Z + F z Uy), 



by Z ^Z ) 1 

2(F X U X + FyUy + F z u z ), 

F x (3u 2 x + u 2 y + u 2 z ) + 2u x (F yU y + F z u z 
F y (u 2 x + 3u 2 y + u 2 z ) + 2u y (F x u x + F z u z 

.2 , „.2 , o„,2x 



F z [u x + u z y + 3< 



2U Z i y F x U X ~\- FyUy J 



F x tiyti z ~\~ Fyti x ti z -\- F z u x u 



y 
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F x u x (3{u 2 y + u 2 z )-A)+FyUy (3(u 2 



+F z u z (3{u 2 x + u\ 
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13F T 



13F, 



y 



13F, 



(G.5) 



Using m%, the source terms in velocity space can be obtained by a procedure 
involving exact inversion that invokes orthogonal properties of the collision 



matrix (see the discussion following Eq. (33) for the D3Q27 lattice). The 
results are summarized as follows: 



So — 




5mg + m 8 4 ] , 
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+ 45m?, + 15mg 


- I0m s 9 


-9m 8 


-m 8 14 
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- 18m 8 
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-m 8 14 
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is |12fflS 
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- 30m 8 - 10m 8 , 


- 9m 8 12 


- m 8 14 ] , 





48 



s 7 = 



— [12m* 
180 L 

' [48m* + 72m* + 72m*. + 

90m* 3 + m* 4 ] , 



720 



18m 3 - 30m*, - 10m* + 9m* 2 
72m* + 90m| 

■ s n + 9m* 2 + 
+ 72m; 

-9m* + 9771^ + 9 ^i2 _ 90mi 3 + m° u \ 



m 



14J > 



+ 90m*- + 90m" + 40m' 



90m* + 90m* 



+ 40m* 



90m* - 90m* 



+ 40m*. 



+9mJ + 9m s n 
S 8 = -^- [48m* - 72m* + 72m* + 72m* - 90m 4 

- 9m* x + 9m^ 2 — »u/« 13 
S g = — [48m* + 72m* - 72m* + 72m* - 90ra| + 90m* - 90m s 6 + 40m 

+9mi - 9m s n + 9m* 2 - 90m* 3 + m s u ] , 
S w = -L [48m* - 72m* - 72m* + 72m* + 90m s 4 

-9m* - 9m s u + 9m* 2 + 90m^ 3 + m* 4 ] , 
S n = — [48m* + 72m* + 72m*. - 72m* + 90m| 
+9mi + 9m*! - 9m* 2 - 90m* 3 + m s u ] , 

5 12 = -L [48m* - 72m* + 72m* - 72m* - 90m* 

- 9m* 2 + 90m* 3 + m* 4 ] , 

5 13 = [48m* + 72m* - 72m* - 72m* - 90m* - 90m* + 90m^ + 40m^ 
+9mi - 9m s n - 9m* 2 + 90m* 3 + m s u ] , 

- 90m s A + 90m* + 90m; 



720 

-9m* + 9m\ 
- [48m*, 



- 90m* - 90mg + 40m* 



+ 90m* - 90m* + 40m* 



- 9m* 2 + 
1 - [48m* - 72m* 
-9m* 



:* - 72m* - 72m* + 1 



• S|! 720 ' 

- 9m s u - 9m* 2 - 90mi 3 + m* 4 ] . 



+ 40m*, 
(G.6) 



For obtaining explicit expressions for the collision kernel, it is convenient to 
express the non-conserved transformed raw moments using the operator no- 



tation given in Eq. (31), which are given as subsets of the particle velocity 



directions for the D3Q15 lattice. It follows that 



K xy ~ {fa\ e ax€ 



ay I 



K yz ~ {fa\ e ay&c 



14 

/ a e ax^ay = 

a=0 


( M 

(? 


-E 
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14 


:? 


B 5 

-E 
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f a^ay^az = \ 

a=0 


( A 6 
V a 


B 6 

-E 

a 



1 fa, 
fa, 
f a) 



(fa\ e ax) — E fa e ax ~ E I ® -^a' 



49 



a=0 \ a / 

14 - / ^ \ 

K zz = (f a \el z ) = /a e L = S ® /a, 
a=0 \ a J 

^ 14 /Aio Bio\ 

K :ct/2/ = (/al e aa: e Qj;) = 5Z f a e ax^ ay = I 5Z _ 5Z I ® /a' 

a=0 \ a a J 

_ 14 _ /An #ii\ _ 

= (f a\ e ax e ay) = 5Z f a e ax e ay = I 5Z ~~ 51 I ® /a' 

a=0 \ a « / 

_ 14 _ /Ai 2 B 12 \ _ 

= (fa\ e ax e az) = 51 f a e ax e az = I 51 _ 51 I ® /a) 

a=0 \ a a / 

_ 14 _ /Ais Bi 3 \ _ 

(/al e aa! e ay e az) = 5Z f a^ax^-ay^az = I 5Z ~~ 5Z I ® /a' 

_ 14 _ / j4l4 \ - 

(/o;l e ax e a2/) = 51 f a e ax e ay = I 5Z I ® /a' (G-7) 

a=0 V a J 



n xxy 



xyz 



n xxyy ' 



where 



A 4 = {7, 10, 11, 14}, S 4 = {8, 9, 12, 13}, 
= {7, 9, 12, 14}, ^ = {8, 10, 11, 13}, 
A 6 = {7, 8, 13, 14}, £ 6 = {9,10, 11, 12}, 
A 7 = {1, 2, 7, 8, 9, 10, 11, 12, 13, 14} , 
A 8 = {3, 4, 7, 8, 9, 10, 11, 12, 13, 14} , 
A 9 = {5, 6, 7, 8, 9, 10, 11, 12, 13, 14} , 
A w = {7, 9, 11, 13} , B w = {8, 10, 12, 14} , 
A u = {7, 8, 11, 12} , B u = {9, 10, 13, 14} , 
A 12 = {7, 8, 9, 10} , B 12 = {11, 12, 13, 14} , 
A 13 = {7, 10, 12, 13}, B 13 = {8, 9, 11, 14}, 
A u = {7, 8, 9, 10, 11, 12, 13, 14}. 



G.3 Collision Kernel 



Following the same procedure and the notations as used for the D3Q27 lat- 
tice and considering factorized attractors for the higher order moments, the 
cascaded form of the central moment collision operator in the presence of 
forcing terms can be constructed. The results are summarized as follows (for 
collisional invariants, g = g\ = g 2 = g 3 = 0): 
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94 
95 
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97 
98 

99 

9w 
9u 
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913 = 

914 = 



U 4 
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u 6 



ixy 



pU X Uy+ -{O x Uy + d y U XJ 



-Vxz + P U xUz + ^ip'x U z + K U x) 



Vyz + P U y U z + ^{VyUz + O z U y ) 
(Ix - %y) + P( U l - U D + (?x U s - ®y^y) 

(ix + ^ - 2 ^J + p(u 2 x + u 2 y - 2u 2 z ) 
+(a x u x + a y u y - 2a z u z ) 



co 7 

T 

12 
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+ (a a ,-u :r + dyUy + + p 
16 



1 , 



~V X yy + 2Uy71 xy + U x T] yy - 2 P U x U y - -B x U y - B y U y U x 



+u y g 4 + -u x (-g 7 + g s + 3<? 9 ), 
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Ixxy 



2u x 7] + u,,,,.,. - 2pu x u y 



-O y U x ~ O x U X Uy 
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-Vxxz + 2u xVxz + UzVa 
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2pu\u z 



1~' 2 



a x u x u z 



+u x g 5 + -u z (g 7 + g 8 + Sg g ), 
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Vxyz + U xVyz + U yV X z + U zVxy 

+a y u x u z + a' z u x Uy) 
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2pu x u y u z - - [a x u y u z 



(G.8) 
(G.9) 
(G.10) 
(G.11) 

(G.12) 
(G.13) 

(G.14) 

(G.15) 

(G.16) 
(G.17) 
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+ 2^7? 



xyy 



— 2— 2 

^UyVxxy ~ U xVyy ~ U yV 



y l xx 



4u x ii y r, 



xy 



+K xx K yy + 3pu x u y + a x u x u y + d y u y u x - 2u x u y g 4 u y )g 7 



+ - 



ul)g% + 



g 9 + 2m^io + 2u y g u ,(G.18) 



where W4, W5, . . . , uiu are relaxation parameters. Note that similar to the D3Q27 
lattice, we have the following relation for the shear viscosity of the fluid 
z/ = c s(^~~ l)' wnere u v = Uj and j = 4, 5, 6, 7, 8. The remaining parameters 
can be adjusted independently to control numerical stability. 
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G.4 Operational Steps 



Finally, by expanding the elements of the matrix multiplication of /C with 



in Eq. (16), the post-collision values of the distribution function augmented 



by source terms corresponding to the D3Q15 lattice are 



7o = 7o + [do - 2(?9 + 32<7i 4 ] + Sb, 
fx = 7i + [do + 9i + 07 + 9s - 09 - 

J 2=1 2+ [do ~ 01 + 07 + 08 ~ 

7 3 = 7 3 + [do + 92 - 97 + 98 - 99 - Bg n - 8g u ] + S 3 , 
[do -92- 



3<?io - 8pi 4 ] + Si, 
- <?9 + 8^10 - 8<? 14 ] + S 2 , 



02 - 07 + #8 - 99 + 8^11 - 80 M ] + S 4 , 
7 5 = 7 5 + [<7o + ^3 - 2^8 ~ 09 ~ 8^12 - 8^14] + S 5 , 

'8 - ^9 + 8^12 - 8^14] + S 6 , 
7 7 = 7 7 + [do + 91 +92 + 93 +94 + 95+96 + 



/ 4 =/ 4 

7 5 = 7 

/e = ^6 



H013 + 

7 8 = 7 8 + [<?o 

"013 + 2^m] + S t 



'14] + s 7 , 

0^ + 02 + 03 " 04 " 05 + 



/q = /o + [fo + 0i - £2 + £3 - 04 

-013 + 2^14] + S 9 , 

7io = 7io + [<?o -9i- 



+913 + 

fll = fll 

-913 + 



9x "02 + 03 + 04 

14] + S10, 

[do + 01+02-03 + 04 

'14} + Sn, 



/12 — /12 + [do — 



9i 



^9i3 + 2gu] + S 12 , 



fl3 — fx3 + [do 

+di3 + z> 
7 14 = 7 14 + [do 

-913 + 2; 



d2-g3-g4 + db- 
'12, 

f 01 - 02 - 03 - 04 

h4\ + S13, 

- di - d2 - d3 + 

114] + S14. 



g g + 2g 10 + 2g n + 1 
g 6 + d9- 2<?io + 2011 + i 

dh-d& + d9 + 2^10 - 2^11 + i 
-ds-d6 + d9- 2gio - 20ii + 

- #5 - #6 + #9 + 2^io + 20n - 

g 6 + d9- 2<?io + ~ 
#5 + da + d9 + 2^io - 
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